
% Classical Minimum Distance (CMD) objective function

function f = cmd_objective4(theta,auxdata)

    % Extract vector of sample moments
    mhat        = auxdata{1};
    
    % Extract vector of sample moments from SCF
    mhat_scf    = auxdata{2};
    
    % Extract the fixed income asset share of group 1
    sa1         = auxdata{3};

    % Extract the invV matrix
    invV         = auxdata{4};
    
    % Extract calibrated parameters from scf
    mhat_scf_param4         = auxdata{5};  
    
    % Extract individual parameters from parameter vector
    pi_i1       = theta(1);
    pi_i2       = theta(2);
    pi_c1       = theta(3);
    pi_c2       = theta(4);
    
    % Extract individual sample moments from SCF
    mu_r1       = mhat_scf_param4(1);
    mu_r2       = mhat_scf_param4(2);
    mu_a1       = mhat_scf_param4(3);
    mu_a2       = mhat_scf_param4(4);
    var_r1      = mhat_scf_param4(5);
    var_r2      = mhat_scf_param4(6);
    var_a1      = mhat_scf_param4(7);
    var_a2      = mhat_scf_param4(8);
    
    % Extract individual sample moments from SCF
    cov_r1r2    = mhat_scf(1);
    cov_r1a1    = mhat_scf(2);
    cov_r1a2    = mhat_scf(3);
    cov_r2a1    = mhat_scf(4);
    cov_r2a2    = mhat_scf(5);
    cov_a1a2    = mhat_scf(6);

    % Model moments
    mu_y1       = mu_r1+mu_a1;
    mu_y2       = mu_r2+mu_a2;
    mu_a        = sa1*mu_a1+(1-sa1)*mu_a2;
    mu_ri       = pi_i1*mu_r1+pi_i2*mu_r2;
    mu_rc       = pi_c1*mu_r1+pi_c2*mu_r2;
    var_y1      = var_r1+var_a1+2*cov_r1a1;
    var_y2      = var_r2+var_a2+2*cov_r2a2;
    var_a       = (sa1^2)*var_a1+((1-sa1)^2)*var_a2+2*sa1*(1-sa1)*cov_a1a2;
    var_ri      = (pi_i1^2)*var_r1+(pi_i2^2)*var_r2+2*pi_i1*pi_i2*cov_r1r2;
    var_rc      = (pi_c1^2)*var_r1+(pi_c2^2)*var_r2+2*pi_c1*pi_c2*cov_r1r2;
    cov_y1y2    = cov_r1r2+cov_r1a2+cov_r2a1+cov_a1a2;
    cov_y1a     = sa1*cov_r1a1+sa1*var_a1+...
                (1-sa1)*cov_r1a2+(1-sa1)*cov_a1a2;
    cov_y1ri    = pi_i1*var_r1+pi_i1*cov_r1a1+...
                pi_i2*cov_r1r2+pi_i2*cov_r2a1;
    cov_y1rc    = pi_c1*var_r1+pi_c1*cov_r1a1+...
                pi_c2*cov_r1r2+pi_c2*cov_r2a1;
    cov_y2a     = sa1*cov_r2a1+sa1*cov_a1a2+...
                (1-sa1)*cov_r2a2+(1-sa1)*var_a2;
    cov_y2ri    = pi_i1*cov_r1r2+pi_i1*cov_r1a2+...
                pi_i2*var_r2+pi_i2*cov_r2a2;
    cov_y2rc    = pi_c1*cov_r1r2+pi_c1*cov_r1a2+...
                pi_c2*var_r2+pi_c2*cov_r2a2;
    cov_ari     = pi_i1*sa1*cov_r1a1+pi_i1*(1-sa1)*cov_r1a2+...
                pi_i2*sa1*cov_r2a1+pi_i2*(1-sa1)*cov_r2a2;
    cov_arc     = pi_c1*sa1*cov_r1a1+pi_c1*(1-sa1)*cov_r1a2+...
                pi_c2*sa1*cov_r2a1+pi_c2*(1-sa1)*cov_r2a2;
    cov_rirc    = pi_i1*pi_c1*var_r1+(pi_i1*pi_c2+pi_i2*pi_c1)*cov_r1r2+...
                pi_i2*pi_c2*var_r2;

    % Vector of model moments
    m           = [mu_y1;mu_y2;mu_a;mu_ri;mu_rc;var_y1;cov_y1y2;...
                cov_y1a;cov_y1ri;cov_y1rc;var_y2;cov_y2a;cov_y2ri;...
                cov_y2rc;var_a;cov_ari;cov_arc;var_ri;cov_rirc;var_rc];
            
    % Difference between model and sample moments
    d           = m-mhat;
    
    % Objective function
    f           = d'*invV*d;
    